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Abstract 

Shear deformation of partially molten rock in laboratory experiments causes the emer¬ 
gence of melt-enriched sheets (bands in cross-section) that are aligned at about 15-20° 
to the shear plane. Deformation and deviatoric stress also cause the coherent alignment 
of pores at the grain scale. This leads to a melt-preferred orientation that may, in turn, 
give rise to an anisotropic permeability. Here we develop a simple, general model of 
anisotropic permeability in partially molten rocks. We use linearised analysis and nonlin¬ 
ear numerical solutions to investigate its behaviour under simple-shear deformation. In 
particular, we consider implications of the model for the emergence and angle of melt-rich 
bands. Anisotropic permeability affects the angle of bands and, in a certain parameter 
regime, it can give rise to low angles consistent with experiments. However, the condi¬ 
tions required for this regime have a narrow range and seem unlikely to be entirely met 
by experiments. Anisotropic permeability may nonetheless affect melt transport and the 
behaviour of partially molten rocks in Earth’s mantle. 

1 Introduction 

The dynamics of partially molten regions of the mantle are not well nnderstood; in part, 
this is dne to the inaccessibility of these regions to direct observation. It is widely agreed, 
however, that laboratory experiments performed on synthetic, partially molten aggregates of 
olivine and basalt [e.g Holtzman et ah, 2003, King et ah, 2010] are an effective test of theoretical 
models of magma dynamics [e.g. McKenzie, 1984]. One featnre of the experiments that remains 
challenging to reprodnce in models is the emergence of high-porosity, melt-rich bands at angles 
between 15° and 20° to the shear plane [King et ah, 2010]. Althongh some characteristics of 
the band pattern are sensitive to the applied strain rate [Holtzman and Kohlstedt, 2007], the 
low angles of these bands are robnst to variation in experimental parameters. Here we consider 
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whether the low angles could arise as a consequence of a permeability that is directionally 
dependent (i.e. anisotropic). 

Anisotropic permeability could be a consequence of another empirically known feature of 
partially molten rocks subjected to deviatoric stress: the microstructural alignment of intercon¬ 
nected pockets of melt between solid grains. This is called melt-preferred orientation (MPO) 
and has been observed in many laboratory studies [e.g. Bussod and Christie, 1991, Daines and 
Kohlstedt, 1997, Takei, 2010]. The alignment may be attributed to the instantaneous state of 
deviatoric stress [Daines and Kohlstedt, 1997, Takei and Holtzman, 2009a], or to the combined 
effects of finite strain, lattice-preferred orientation, and anisotropic surface energy of olivine 
grains [Bussod and Christie, 1991, Daines and Kohlstedt, 1997, Jung and Waff, 1998]; it is 
likely some combination of the two. Since the Darcian permeability of partially molten rocks 
arises from the shape and interconnectedness of melt pockets at the grain scale [e.g. Bear, 
1972, Scheidegger, 1974], it is reasonable to assume that the anisotropic alignment of pores 
between grains leads to anisotropy in permeability. Daines and Kohlstedt [1997] estimated 
this anisotropy as a function of differential stress and found that permeability in the direction 
parallel to the maximum compressive stress cxi was enhanced by a factor of up to 15 over that 
parallel to the direction of maximum tensile stress. This is consistent with a theoretical model 
for anisotropy of permeability due to MPO by Hier-Majumder [2011]. 

Since both melt-banding at the macroscopic scale and melt-preferred orientation at the 
microscopic scale emerge under the same physical conditions, it is logical to ask whether their 
dynamics are linked. In particular, the question we address here is whether the low angle 
of high-porosity bands observed in experiments [see schematic diagram in Takei and Katz, 
2013, fig. 1] can be explained through consideration of anisotropic permeability arising from 
MPO. We therefore develop and analyse two-phase models based on viscous compaction theory 
[McKenzie, 1984]. Stevenson [1989] showed that for a matrix viscosity that weakens with 
porosity, unstable growth of porosity perturbations can lead to localisation of melt. We adopt 
this rheology and modify the two-phase theory with assumptions of how deviatoric stress might 
modify the permeability of the aggregate. Calculations show that anisotropic permeability does 
indeed exert control over the predicted angle of melt-rich bands; it can even give rise to the 
low angles observed in experiments. However, the parametric conditions required to reproduce 
the angles in experiments are a rather restrictive set, making it unlikely that this is the true 
explanation of the observations. 

Other theoretical approaches have been made to explain the emergence and angle of melt- 
rich bands. In all cases besides the present one, authors have introduced more complicated 
rheological formulations to obtain low band angles. Katz et ah [2006] found that low angles 
are predicted under non-Newtonian viscosity, but experiments by King et ah [2010] produced 
low-angle bands and had a measured stress exponent of ~1. A more recent approach considers 
the effect of deviatoric stress and MPO on diffusion-creep (Newtonian) viscosity [Takei and 
Holtzman, 2009a,b]. Melt at the grain boundaries is a fast pathway for diffusion of the com¬ 
ponent that makes up solid grain; coherently aligned melt pockets thus give rise to anisotropy 
of aggregate viscosity under diffusion creep [Takei and Holtzman, 2009a] . The consequences of 
microstructural and viscous anisotropy for melt-band formation were investigated by Takei and 
Holtzman [2009c], Butler [2012], Takei and Katz [2013], Katz and Takei [2013], Allwright and 
Katz [2014], and Takei and Katz [2015]. All of these studies predict low-angle, high-porosity 
bands, consistent with experimental results, but none of them address the implications of MPO 
for permeability and melt segregation. 

The manuscript is organised as follows. In §2 we present the governing conservation equa¬ 
tions, the constitutive laws and a useful rescaling. §3 presents the assumptions and formulation 
of the tensorial permeability. In §4 we develop the linearised stability analysis. This section 
contains sub-sections considering the effect of the direction of MPO (§4.1), wavelength of poros- 
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ity perturbations (§4.2), and the anisotropy conditions that can give rise to low-angle bands 
(§4.3). Results from numerical models are presented in §5 and compared with results from the 
stability analysis. In §6 we discuss, summarise, and conclude the manuscript. 

2 Governing equations 

2.1 Conservation statements 

Equations describing the dynamics of an aggregate of solid mantle and liquid magma have been 
derived by various authors [e.g. McKenzie, 1984, Bercovici et ah, 2001, Rudge et ah, 2011]. Here 
we consider a simplified form of the equations that is relevant to laboratory experiments. In 
this case there is no melting and deformation is rapid enough that we can neglect the effect of 
body forces. The equations are written in terms of liquid-volume fraction 0, liquid velocity v^, 
liquid pressure p^, and matrix velocity as 


dt(t) + V ■ [0u^] , 

(la) 

001 - 0) + V ■ [(1 - 0)uT , 

(lb) 

0 {v^ - v^) + (K//i) Vp^, 

(Ic) 


(Id) 


where p is the liquid viscosity, W = -|- (1 — (j))cr^ is the stress tensor of the two-phase 

aggregate, and K is the permeability tensor. These equations represent conservation of mass 
for the liquid and solid, and balance of forces for the liquid and the aggregate, respectively. 
They incorporate the assumption of constant and uniform phase densities. 

2.2 Constitutive relations 

To close the system of partial differential equations (1), we require the specihcation of consti¬ 
tutive relations for stress and permeability. 

The total stress can be written in terms of the pressure and a deviatoric stress tensor, 
cr^ = —p^l and cr^ = —p^l + , where we have made the usual assumption that ~ 0 . 

Following Rudge et al. [2011], we express the constitutive relationships for compaction and 
shear as 


= (2a) 

1-0 

[Vn^ + (Vn^)^ - fiV ■ n^] . (2b) 

1 — 0 ^ ^ 

This is consistent with the formulation of McKenzie [1984] and gives 

a = -p^l + [Vn^ + (Vn^)'^ - fiV ■ v^] + IQV ■ (3) 

as the phase-averaged stress tensor, and are the aggregate shear and compaction viscosity. 
For present purposes, it is sufficient to take forms for these that are theoretically justihed 
[Simpson et ah, 2010a,b] and minimally complicated [Stevenson, 1989]: p^ = exp[—A(0 —0o)] 
and C(/> = the shear viscosity, A controls the reduction in viscosity due to porosity; in 

the bulk viscosity, controls the ratio of compaction to shear viscosity. Parameters Pq and 0o 
are constant, reference values. Stevenson [1989] showed that the porosity-weakening of viscosity 
(A > 0) gives rise to a melt-banding instability. Here we take A = 27, = 5/3, and 0o = 0.05 

[e.g. Kelemen et ah, 1997, Takei and Holtzman, 2009a, Holtzman et ah, 2003, respectively]. 


3 




The permeability can be broken down into an isotropic factor, dependent on the porosity, 
and an anisotropic factor that we use to model the effect of melt-preferred orientation. This 
gives 

K = Ko(0/0o)'A, (4) 

where Kq is a reference permeability, £ is a constant measured to be between two and three [e.g. 
von Bargen and Waff, 1986, Miller et ah, 2014], and A is a non-dimensional, second rank tensor 
determined by the form of the permeability anisotropy. Note that A could depend on other 
variables, particularly the deviatoric stress or hnite strain that cause MPO. Here we neglect 
these dependencies and take £ = 3 for consistency with previous studies. 

2.3 Rescaling 

Equations (1), (3), and (4) can be combined to eliminate resulting in the system 

dt(j) = V ■ [(1 - 0)v], (5a) 

V • V = V ■ [(K//i)VP], (5b) 

VP = V ■ [r 7 <^(Vv + Vv^)] + V [(C^ - V ■ v] , (5c) 

where we have substituted P p^, \ ^ for simplicity of expression. 

We now introduce the following characteristic scales: [a;] 5, [f] ~ 7 ^ [ p^] ^ 7 o, [v] 7 ( 5 , 

[P] ~ 707 , and [K] Kq. Here we have used the compaction length 5 = (r^ -|- A/3)poKo/p, 

a characteristic lengthscale for liquid/solid interaction [McKenzie, 1984]. Our focus will be 
simple-shear flows with initially uniform strain rate 7 , which gives us a characteristic time- 
scale for the problem. Using these scales to nondimensionalise all symbols we obtain 


dt(f) = {I - (f))C■ V(j), ( 6 a) 

C = (rc + 4/3)-iV • [(0/0o)'AVP] , ( 6 b) 

VP = V ■ [^^(Vv Vv"^)] -F V [p^ (r^ - 2/3) C ], ( 6 c) 

where we have dehned the compaction rate as C = V ■ v. 

3 Formulation of anisotropic permeability 

The dynamic origins of melt-preferred orientation remain poorly understood. In laboratory 
deformation experiments, melt pockets become elongated such that in cross section, they have 
a long axis and a short axis, as shown in Fig. 1. Coherent alignment of the long axes of 
melt pockets results from an unknown combination of the instantaneous stress [Takei and 
Holtzman, 2009a] and the hnite strain [Daines and Kohlstedt, 1997, Jung and Waff, 1998]. The 
latter acts to align olivine grains, which have anisotropic surface energy and hence anisotropic 
wetting affinity. If the deviatoric stress is the dominant forcing of MPO, we would expect melt 
pockets to be in instantaneous alignment, with their long axes perpendicular to the direction 
of maximum tension (the Us-direction). We can generalise this, however, and investigate the 
consequences of MPO alignment in terms of the angle 0 between the shear plane and the 
normal to the long axis of a melt pocket. 

This angle is used to dehne a rotated coordinate system as in Fig. 1. The a:'-direction is 
normal to melt pockets and expected to have a low permeability. The ^/'-direction is parallel 
to melt pockets (and in the plane dehned by the eigenvectors corresponding to cxi and as); this 
direction is expected to have a high permeability. 
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Figure 1: Schematic diagram of melt-rich bands and melt-preferred orientation, (a) shows a 
macroscopic image of melt-rich bands, aligned at an angle 9 to the shear plane. The image is 
adapted from Figure 4 of King et al. [2010]. A schematic illustration of the simple-shear flow field 
is shown at left, (b) shows a schematic of melt-preferred orientation aligned at an angle 0 to the 
shear plane. The image is adapted from Figure 2 of Holtzman et al. [2003] by thresholding the 
greyscale such that melt is grey and solid is white. The melt-preferred orientation is aligned with a 
rotated coordinate system (x',y') and has enhanced permeability in the direction of the long axes 
of the melt pockets (y'). The scale of this image is microscopic — of order 10 x the grain size — so 
well below the continuum scale of high-porosity bands. 


We quantify these expectations with the tensor A' in the rotated coordinate system, 


A' 


k± 0 
0 h 


(7) 


where > 1 is the permeability enhancement along y' and 0 < < 1 is the reduction along 

x'. A symmetry argument (Appendix A) shows that A' must be diagonal. Then, by rotating 
from the primed to the unprimed coordinates, the anisotropy matrix is written as 


A = RA'R^ 


/cj_ cos^ 0 -|- k\\ sin^ 0 (/c_i_ — k\\) cos 0 sin 0 
(fc_i_ — k \\) cos 0 sin 0 fc_i_ sin^ 0 -|- k\\ cos^ 0 


( 8 ) 


Here we could factor out fey, for example, and lump it with Kq, reducing by one the number of 
parameters in the problem. However, since we have non-dimensionalised lengths with the com¬ 
paction length (and hence with Kq), this would obscure the role of permeability enhancement 
parallel to MPO. 

Much of the work below is in understanding the behaviour of solutions to the system (6) 
under different assumptions for the values of kj_, k\\, and 0 in the tensor (8). A list of the key 
non-dimensional symbols used in this manuscript is provided in table 1. 


4 Linearised stability analysis 

Equations (6) form a coupled, non-linear, time-dependent system. To make analytical progress, 
we follow previous workers [e.g Spiegelman, 2003, Katz et ah, 2006] and employ a linearised 
stability analysis. We expand the problem variables in a power-series of e <C 1 and truncate 
after Erst order, 

{ 0 = 00 + {x,t) , 

P = Po + ePi{x,t), 

V = (a;)-|-(a;, t), 

C = Cq T • 

Terms of order e*’ are taken to dehne the base-state about which we perturb. The expansion of 
porosity can be used to express the nondimensional permeability and viscosity to hrst order as 
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Symbol 

meaning 

equation, value, or range 

X 

position 


t 

time 


P 

liqnid pressure 


V 

velocity 


0 

porosity 


C 

compaction rate 

V • V 


shear viscosity 

1 ]^ = exp [-A(0 - 0o)] 

a 

bnlk viscosity 

C<A = rcV4> 

A 

porosity-weakening exponent 

27 


viscosity ratio 

5/3 

00 

reference porosity 

0.05 

1 

permeability exponent 

3 

A 

permeability anisotropy tensor 

eqn. ( 8 ) 

fc|| 

MP 0-parallel permeability factor 

1-200 

/C_L 

MPO-perpendicular permeability factor 

0.005-1 

0 

angle of MPO to shear plane 

0°-90° 

9 

angle of pertnrbation wavevector to shear plane 

0 

0 

00 

1 

0 

0 

e 

very small number 

< 1 

s, s 

perturbation amplitnde, growth rate 

eqn. (12), (13) 

K, K 

pertnrbation wave-vector, magnitnde 

K = K (sin 6 *, cos 6 *) 


Table 1: Dimensionless symbols and their meaning. 


( 0 / 0 o)^ = 1 + ei(j)i/(t)o ( 10 a) 

g-A(<^-0o) = X _ (xOb) 


Substituting (9) and (10) into ( 6 ) and requiring terms to balance at each order allows us to 
derive two systems of linear eqnations at C>(e°) and We assume that the leading-order 

porosity 0o is nniform and that the flow is simple shear with = ( 7 ?/, 0). The uniform, 
leading-order viscosity and permeability therefore resnlt in Co = 0 and Pq = const. 

The first order balance of equations ( 6 ) becomes 

at 0 i = (l- 0 o)Ci-v(°i-V 0 i, ( 11 a) 

Ci = (r^ + 4/3)-'V-[AVPi], (11b) 

VPi = V ■ + (Vv(^i)^] -b (rc - 2/3)VCi - AV ■ [(j)i (Vv(°) -f (Vv^^i)^)] . (11c) 


We consider harmonic porosity perturbations that are advected by the backgronnd simple shear 


01 = exp 



( 12 ) 


where exps(t) is the time-dependent amplitnde and the wave-vector is given as k = {kx, Hy) = 
K (sin 6*, cos 6*). Hence 9 represents the angle between the pertnrbation wave vector n and the 
?/-axis or, eqnivalently, between the perturbation wavefronts and the shear plane. 

This linearised system (11) can be solved, along with (12), to obtain the growth rate (details 
in Appendix B) 


s = 


1-00 
+ 4/3 


1 -b (k, Ak) 



(13) 
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where the inner prodnct notation means {k, Ak) = X] S 

i j 

For isotropic permeability kj_ = k\\ = 1] then from eqnation ( 8 ), A rednces to the identity 
matrix. In this case, the growth rate is identical to that calcnlated by Spiegelman [2003] for 
isotropic permeability (and viscosity). Fnrther, in the limit of vanishing pertnrbation wave¬ 
length K —)■ cx) (where k = |k|), the ratio {k,Ak) /{1 + (k, Ak)) —)■ 1 and so the growth rate 
again rednces to that of the isotropic case. In contrast, the limit k = 0 has zero growth rate. 
It can be shown (Appendix B) that the greatest sensitivity to the anisotropy of permeability 
occurs where 271 /k ~ 1 or, in dimensional terms, when the wavelength of the instability is 
approximately equal to the compaction length. 



Figure 2: Anisotropy factor F as a function of 0-1-0 for k = 1 and various values of fey. (a) k± = 0.5. 
(b) k± = 0.005. 


For a permeability matrix A from equation (8) we can expand and simplify to obtain 

(k, Ak) = [/c|| — (^11 — k±) sin^(0 -|- 0)] . 

Hence from equation (13), we are interested in the behaviour of the factor 

^ ^ h - (fc|| -fcu)sin ^(0 + 0 )] 

1 + K? [k\i — {k\\ — k±) sin ^(0 -|- 0 )] 

that modifies the isotropic growth rate. F is plotted in Fig. 2 as a function of 0 -|- 0 for k = 1 
and various values of k\\ and MPO is parallel to porosity perturbations when 6 + Q = 90° 
or 270°. At these two angles, the reduced permeability kj_ in the MPO-perpendicular direction 
controls the flow of melt perpendicular to bands and causes a minimum in F. With k± —?• 0, 
melt cannot segregate in the direction perpendicular to MPO, and hence it cannot feed bands 
that are oriented parallel to MPO. When 9 + Q ^ 90° or 270°, values of /cy > 1 facilitate 
segregation of melt into (or out of) perturbations. 

Although permeability can promote or impede melt transport, it does not drive melt trans¬ 
port. In this context, it is important to note that in the linearised analysis, the permeability 
does not multiply any zeroth-order fields and hence perturbations to the permeability play no 
role in the first-order balance (11). Only anisotropy of the zeroth-order permeability affects the 
growth rate of porosity perturbations. Here we assume that MPO (and hence anisotropic per¬ 
meability) emerges on a much shorter time scale than melt redistribution, so the leading-order 
MPO is in equilibrium with the leading-order stress state. This assumption is consistent with 
grain-scale models [Hier-Majumder, 2011] as well as experimental observations [Takei, 2010, 
Zimmerman et ah, 1999]. Other experiments indicate that the strength of MPO increases with 
the amount of time that the stress is applied [Daines and Kohlstedt, 1997]. Furthermore, the 
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local stress state is affected by the segregation of melt into bands [e.g. Takei and Katz, 2015], 
which could modify the MPO and the permeability. However, in the context of linearised 
analysis, this would not appear at leading order and hence would not modify our results. 

Equation (13) shows that the growth-rate of perturbations depends on the product of F 
with sin26', where the latter represents the deviatoric stress (tension positive) that drives melt 
segregation. In the absence of anisotropic permeability, sin26' is the sole angular dependence 
of the growth rate and bands are predicted to grow fastest at 45° to the shear plane. We next 
explore how and when F modihes this prediction. 

4.1 MPO orientation and growth rate 

We return to (13) and consider the role of the orientation of anisotropy 0. Within the lim¬ 
itations of linearised analysis, it is the porosity perturbation with the largest instantaneous 
growth rate that is predicted to dominate after hnite time, and hence to be expressed by the 
physical system. We are therefore interested in the angle 9* at which the growth rate reaches 
its maximum, s*. Since the growth rate is most sensitive to anisotropy at k ~ 1, we consider 
only K = 1 in this section; this constraint is relaxed in the following sections. 



Figure 3: Growth-rate s versus 6 for various values of anisotropy angle 0. Each curve is calculated 
with k± = 0.005, fc|| = 2. 

Fig. 3 shows the growth rate versus 6 for a set of values of anisotropy angle 0. Indepen¬ 
dent of any imposed anisotropy, perturbation growth is driven by tension across a high-porosity, 
viscously weak band [Stevenson, 1989]. As band angles rotate into directions parallel or perpen¬ 
dicular to the shear plane, this tension vanishes and the growth rate goes to zero [Spiegelman, 
2003]; this is why all curves in Fig. 3 cross s = 0 at 0° and 90°. Anisotropy of permeability mod¬ 
ihes the sin 20 dependence of growth-rate predicted by Spiegelman [2003]. As 0 is increased 
from its smallest value, a local minimum in the positive growth rates appears in the range 
0° < 0 < 90°. This minimum exists at the angle where perturbation wavefronts are parallel to 
MPO, 0 = 90° — 0. For the case of 0 = 45°, the minimum occurs at 0 = 45° and we hnd two 
peaks of equal growth rate, one at an angle lower than 45° and one higher. Katz et al. [2006] 
show that although the two perturbation angles have the same instantaneous growth rate, it is 
the lower-angle perturbation that dominates, since the high-angle band is rotated more rapidly 
by the background, simple-shear flow out of favourable orientation for growth (eqn. (12)). 








4.2 Perturbation wavelength and growth rate 

As in previous work on isotropic permeability [e.g. Stevenson, 1989], the growth rate in equa¬ 
tion (13) depends on the wavelength of perturbations relative to the compaction length. Fig. 4(a) 
plots growth rate against the log of perturbation wavelength for the angle 6* that has maximum 
growth rate. Here, we consider only MPO aligned at 0 = 45°. For perturbation wavelengths 
much larger than the compaction length, the growth rate of perturbations drops to zero. Over 
these distances, pressure gradients associated with variation in viscosity are small and cannot 
drive melt segregation and differential compaction. 



Figure 4: The effect of perturbation wavelength 27r/K when 0 = 45°, k± = 0.005 and fey = 2. 
(a) Maximum growth rate, s*, against log-wavelength. (b) 9* against log-wavelength. The angle of 
both the low- and the high-angle growth peaks is plotted. 

Looking instead at wavelengths much smaller than the compaction length. Fig. 4(a) indicates 
that there is no wavelength selection by the growth rate. At vanishingly small perturbation 
wavelengths, the resistance to segregation by Darcy drag is negligible. This highly local seg¬ 
regation is insensitive to the permeability; its orientation of maximal growth rate is therefore 
determined entirely by the orientation of the stress tensor and hence by sin 20. Fig. 4(b) plots 
6* as a function of perturbation wavelength; it confirms that at wavelengths smaller than a 
tenth of the compaction length 6* = 45°, as in the isotropic case. For increasing wavelength, 
anisotropy of permeability plays an increasing role; the 45° growth-rate peaks splits into two 
peaks of equal height (as in Fig. 3 for 0 = 45°). However, this split occurs as the growth 
rate declines to zero. Therefore, in a full, nonlinear solution to the governing equations that 
is initialised with a broadband spectrum of perturbations, the growth of shorter-wavelength 
perturbations would dominate the porosity held. These would likely appear at ~ 45°. However 
Butler [2010] showed that the wavelength of high porosity bands in numerical simulations in¬ 
creases as the system evolves, due to the kinematics of the background how. He explained that 
bands grow, rotate, and widen until wavelengths are similar to one compaction length. This is, 
as we have shown, the length scale at which anisotropic permeability has a signihcant effect on 
band angles. 
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It is also important to note, however, that application of the present model to arbitrarily 
small wavelength is invalid because the continuum assumptions hold only at length-scales much 
larger than the solid grains. Indeed if we consider the formation of “bands” at the grain scale, 
the distribution of melt is indistinguishable from melt-preferred orientation. So for a variation 
in porosity to be considered a melt-rich band, the variation must take place over a length-scale 
on which the concept of porosity is well defined; i.e. a large multiple of the mean grain size. 
Such wavelength are, in fact, observed in experiments [Holtzman et ah, 2003, King et ah, 2010] 
indicating that either the initial variations of porosity are concentrated at intermediate or larger 
scales, or that the non-linear system is regularised over finite time by some process that has not 
been included in the theory above. Possible examples of the latter are surface-energy-driven 
flow [Bercovici and Rudge, 2015] and dissolution/transport/precipitation due to gradients in 
chemical potential [Takei and Hier-Majumder, 2009]. 

4.3 The conditions for low-angle bands 

To search for anisotropy conditions conducive to formation of low-angle bands, we map out 
the effect of anisotropy on 9*. Fig. 5 displays contours of 9* and is coloured according to 
the corresponding (maximum) growth rate s*, for fixed k = 1. This is done in the space of 
(— logfc_L, logfc||) because k± ranges between zero and one, while fey is greater than or equal to 
one. For any point in this space at which s has two maxima with equal rates, we select the one 
at smaller 9, corresponding to bands that are more slowly rotated by the base-state flow. In 
these plots, the lower-left corner corresponds to the isotropic case (where 9* = 45°) while the 
top-right corner corresponds to the case with strong anisotropy. 



Figure 5: Contours of 9* (degrees) superimposed on a color-scale plot of s* for k = 1. (a) 0 = 20°. 
(b) 0 = 45°. (c) 0 = 70°. 

Fig. 5 demonstrates that 9* is most effectively reduced by a sharp decrease of k± with little 
or no change in k\\. Panel (b) uses 0 = 45° and shows that in this case, 9* can be lowered 
to below 20° by imposing k ~ 1, fey ~ 1, and k± —)■ 0. Under these conditions the predicted 
angles are consistent with bands observed in experiments, however s is low and so rotation due 
to the background flow is likely to result in bands observed at higher angles than the angle of 
maximal instantaneous growth rate. Panels (a) and (c), however, show that for 0 = 20° and 
0 = 70°, 9* cannot be lowered to angles consistent with the experimental band angles. 

All three panels of Fig. 5 show that away from the origin, the maximal growth rate s* 
increases with fcy and is nearly independent of k±. Close to the origin, where the permeability 
is approximately isotropic, the growth rate does depend on k± (see especially panel (b)). For 
bands to grow, the melt segregation velocity must have a component normal to the pertur¬ 
bation wavefronts. And as we have seen above, growth rate is sensitive to the orientation of 
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perturbations with respect to MPO and to the perturbation wavelength relative to the com¬ 
paction length. A simple heuristic that combines these dependencies is the non-dimensional 
compaction length in the direction normal to perturbation wavefronts of the fastest-growing 
perturbation. We define this as 

5* = \J (k*, Ak*) = K\Jk\\ — {k\\ — k±) sin^( 0 * + 0 ). 

This dehnition requires that when MPO is parallel to maximally growing perturbations (0* -|- 
0 = 90° or 270°), the compaction length across bands is controlled by a/^, whereas when 
MPO is normal to these perturbations, the compaction length is controlled by Returning 

to Fig. 5, far from the origin we have k\\ 3> k± and sin^(0* -|- 0) 7 ^ 1, so 5* (and hence 
s*) is primarily dependent on k\\. Near the origin, k\\ /c_i_ 1 and so 5* depends on both 
permeability parameters. In panel (b) near the origin, we have the special case of 6 ** = 0 = 45°; 
hence S* = \/k± so s* depends only on k± here. 


5 Numerical simulations 

The linearised stability analysis presented above is valid at infinitessimal strains and only for the 
plane-wave perturbations considered. Numerical solutions discussed in this section allow us to 
verify and extend the linearised analysis. We solve the system (6) with permeability tensor (8) 
and shear and compaction viscosity as specihed above in §2. The domain is rectangular with 
height H and boundary conditions of tangengial velocity v = (dimensional; positive 

velocity on the top boundary) and impermeability ^ ■ K = 0 imposed on the top and bottom 
boundaries (the domain is periodic in the x direction). Spatial derivatives are discretised on 
a Cartesian, fully staggered grid with square cells. A semi-imp licit, Crank-Nicolson scheme is 
used to discretise time, and the hyperbolic equation for porosity evolution is solved separately 
from the elliptic system in a Picard loop with two iterations at each time-step. The solutions are 
obtained in the context of the Portable, Extensible Toolkit for Scientific Computation [PETSc, 
Balay et ah, 2001, 2004, Katz et ah, 2007]. Full details and references are provided by Katz 
and Takei [2013]. 

Numerical simulations of band formation with isotropic permeability are documented in the 
literature [e.g. Katz et ah, 2006, Katz and Takei, 2013] and not repeated here. Instead, we 
consider simulations in which the anisotropy is hxed at 0 = 45°, k± = 0.005, and k\\ = 2; we 
vary the ratio of the smallest perturbation wavelength to the compaction length. All simulations 
are initialised with the identical porosity held, shown in Fig. 6(a). This held is constructed 
as (f){x,t = 0) = 00 -|- e(f)i{x), where 0o = 0.05 and e = 0.005. (f>i{x) is a smooth, random 
held with unit amplitude, generated by hltering grid-scale white noise to remove variation at 
wavelengths below 77/10. Hence the smallest perturbation wavelength is 27r /= H/{106), 
in non-dimensional terms that correspond to the abscissa in Fig. 4. A suite of simulations at 
diherent values of compaction length 6 was run to compare the distribution of angles for a 
broadband perturbation. 

If the prediction from the linearised stability analysis holds, we expect simulations to pro¬ 
duce bands at low (~20°) and high (~70°) angles when the minimum perturbation wavelength 
is greater than the compaction length, 77/(105) > 1. In contrast, when the minimum pertur¬ 
bation wavelength is smaller than the compaction length, we expect a transition to 6 > 45°. 
These systematics are discernable in Fig. 4. 

Fig. 6(b) shows the porosity field from a simulation with 6 = 0.0577 (and thus log 27i/Kraax = 
0.3) after a strain of 7 = 0.4. The porosity field, which was initially as shown in panel (a) of 
the same figure, has evolved to form low- and high-angle porosity bands. The angles of these 
features are characterised in the power spectrum shown in panel (c). Here, spectral power from 
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Figure 6 : Results from numerical solutions of (6) with permeability tensor (8) on a 1200x600 grid. 
Parameter values are k± = 0.005, = 2, 0 = 45°. The coordinate axes {x,y) have been rescaled 

with S/H, where 5 is the compaction length and H is the domain height, (a) The initial porosity 
field with mean (j)Q = 0.05 and random perturbation of amplitude e = The perturbation is 

white noise, filtered to remove wavelengths smaller than 1/10*'^ the domain height. This field is 
used to initialise all simulations, (b) The porosity field at strain 7 = 0.4 from a simulation with the 
compaction length chosen as 5 = 0.05H such that log27r/Kmax = 0.3. (c) Spectral power binned by 
perturbation wavefront-angle 9 to the shear plane [after Katz et ah, 2006] for the porosity field at 
7 = 0.4 from simulations with values of 6. Curves are normalised to have unit maximum. 


a two-dimensional Fourier transform is binned according to the angle that wavefronts make 
with the shear plane. The black curve represents the case of perturbations with a minimum 
wavelength that is larger than the compaction length. It is peaked at ~22° to the shear plane, 
with a secondary peak at ~75°. This asymmetry is in contrast with the equal heights of growth- 
rate peaks in Fig. 3 for 0 = 45°. However, the relatively smaller height of the high-angle peak is 
expected because under simple-shear deformation, high-angle features are rotated more rapidly 
than low-angle features; this means that they spend less time in optimal orientation for growth 
[Katz et ah, 2006]. 

Fig. 6 (c) compares the angular spectra of simulations with three different values of log 27r/Kmax- 
As the minimum wavelength of perturbation decreases from longer than to shorter than the com¬ 
paction length, the dominant band orientation rotates to higher angle. For log27r/fi:max ~ —1, 
the low- and high-angle peaks have merged into a single peak at about 45° (slightly higher due 
to rotation). This conhrms the prediction of linearised analysis over hnite strains. 

6 Discussion and conclusion 

The foregoing theory and results extend the equations of two-phase, magma/mantle dynamics 
to include anisotropic permeability. Anisotropic permeability is likely to arise as a consequence 
of melt-preferred orientation, the coherent alignment of melt pockets between solid grains. This 
alignment may be a consequence of hnite strain or of deviatoric stress, or some combination of 
the two. Although we don’t model details of the two-phase microstructure [Hersum, 2008], we 
assume that MPO is a consequence of deviatoric stress, develops instantaneously, and is thus 
active at f = 0 [Takei and Holtzman, 2009a] . We postulate a simple anisotropy tensor with two 
eigenvalues representing the MP0-parallel and MP0-perpendicular permeabilities. An angle 
between the shear plane and the perpendicular to MPO then becomes the third parameter in 
the anisotropy model. 

To understand the behaviour of the system of governing equations and constitutive laws, we 
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perform linearised stability analysis and obtain numerical solutions. We explore the {k\\, k±, 0) 
parameter space of the anisotropy model with a particular focus on 0 = 45°, representing MPO 
aligned perpendicular to the maximum tensile direction of the deviatoric stress tensor. 

For perturbations with wavelengths approximately equal to the compaction length, anisotropic 
permeability modihes the angular dependence of the perturbation growth rate. It suppresses 
the growth of high-porosity bands when those bands are in or near alignment with the MPO 
(high-permeability) direction. This is because melt that segregates into MPO-parallel bands 
must overcome the reduced permeability perpendicular to MPO. For 0 = 45°, MPO is aligned 
normal to the direction of maximum tensile stress, and hence the reduced permeability re¬ 
sists band growth at precisely the angle that would be optimal in the absence of anisotropy. 
In this case, the symmetric peaks in growth rate arise at angles where the product of band- 
perpendicular tension and band-perpendicular permeability is maximised. This angle selection, 
however, depends on the wavelength of perturbation. Wavelengths much longer than the com¬ 
paction length do not generate sufficient pressure gradients to cause melt segregation and band 
growth; wavelengths much shorter than the compaction length are unhindered by the extremely 
short pathways of melt segregation (over which they experience the anisotropy), and hence grow 
optimally at 45°. 

Of particular interest in this exploration is the question of whether anisotropy of perme¬ 
ability can give rise to bands of high porosity that are oriented at low angle (< 20°) to the 
shear plane, as observed in laboratory experiments. Our analysis predicts that it is possible, 
under certain conditions (i.e. 0 ~ 45°, 27i/k ~ 1, k\\/k± ^ 1), to obtain bands oriented at 
the low angles that are consistent with experiments. These are fairly restrictive conditions. 
The requirement that initial porosity perturbations be of a wavelength that is greater than 
about one compaction length seems especially signihcant. Also, the angle of MPO that is typi¬ 
cally measured under relevant experimental conditions is closer to 0 = 30°, as opposed to 45°, 
though there is no accepted explanation for this. 

On the other hand, there are arguments for the conditions above to be met, in which case 
anisotropic permeability might play a role in shaping and orienting the high porosity bands. 
Firstly, bands observed in experiments tend not to have wavelengths of a much smaller scale 
than the compaction length [Holtzman and Kohlstedt, 2007]. This suggest that either the initial 
porosity variations are indeed concentrated at a relatively large scale or, as briefly considered 
above, the non-linear system is regularised by a process that occurs at a scale close to the 
compaction length. This could occur by surface tension acting over a diffuse interface [Bercovici 
and Rudge, 2015] or by diffusion of chemical components [Takei and Hier-Majumder, 2009]. 
Secondly, the argument given by Butler [2010] for the growth of band wavelengths suggests 
that the porosity variations would evolve kinematically to a length scale comparable to the 
compaction length and hence the system would become susceptible to the effects of anisotropic 
permeability. 

Although anisotropic permeability in the form discussed here may not independently explain 
the low angle of porosity bands in laboratory experiments, it may still be relevant in shaping 
the dynamics of those experiments, or of partially molten rocks more generally. Anisotropy 
of permeability may, for example, affect the trajectories of rising magma beneath mid-ocean 
ridges [Morgan, 1987] or melt segregation in magma chambers [Bergantz, 1995]; in industrial 
processes, it has implications for macrosegregation during solidihcation [e.g. Yoo and Viskanta, 
1992]. More work is therefore needed to understand the consequences of anisotropy as well 
as the causes. In this latter category: How can melt-preferred orientation be quantitatively 
related to permeability? And, even more fundamentally, what is the energetic or mechanical 
cause for melt-preferred orientation? 


13 


Acknowledgements The authors thank S. Butler and an anonymous reviewer for insight¬ 
ful comments. The research leading to these results has received funding from the Euro¬ 
pean Research Council (ERC) under the European Unions Seventh Framework Programme 
(FP7/20072013)/ERC grant agreement 279925. JT-W was supported for summer 2014 by as 
Research Experience Placement grant from the Natural Environment Research Council of the 
Research Councils UK. RFK is grateful for the support of the Leverhulme Trust. Numerical 
solutions were computed at the Advanced Research Computing facility of the University of 
Oxford. 


Appendices 


A Diagonality of permeability tensor 


To understand why the permeability tensor A must be diagonal when expressed in a coor¬ 
dinate system aligned with MPO, consider equation (Ic) with uniform porosity 0 = 0o and 
permeability given by equation (4), 

00 = —AVp^. (A.l) 

/i 

We now adopt the primed coordinate system of Figs. 1 and 7 that is aligned with the melt- 
preferred orientation. Making no assumptions about the form of the permeability tensor, equa¬ 
tion (A.l) becomes 


^ f k_i_ fcxi 
P- \ kx2 k^\ 


V/, 


(A.2) 


where q = ~ is the segregation flux. For a unit pressure gradient, = (0,1), 

applied in the p'-direction we have 


q =—{kx^x + kiiy) (A.3) 

P 

with an a:'-component q^i proportional to fexj. 

If we reflect in an axis parallel to the p'-direction, as shown in Fig. 7, the pressure gradient 
is unchanged but the cross-MPO segregation q^/ becomes proportional to —fcxi- However, 
the melt-pocket structure in Fig. 1 is unchanged and so the permeability tensor should not 
change. Hence the segregation obtained from (A.2) is still given by equation (A.3) and we have 
kxi = —kxi = 0, and by a similar argument kx^ = 0. 


B Solving the linearised equations for s 

Because the system of equations that governs the perturbation quantities is linear, we are able 
to relate the other variables to 0i via unknown coefficients, 

vW=v0i, Pi = P0i, (B.l) 

The order-e balance derived from equation (6a) then gives 

^ = (1 - 0o) V ■ V01 - ■ V01. (B.2) 
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Figure 7: Reflection of MPO in y'-axis 


Expanding time and space derivatives using (12) and rearranging, equation (B.2) becomes an 
equation for the growth rate 


S = (1 - 0o) ■ V = (1 - 0o) Cl/01- 


(B.3) 


Substituting = (E, 0), = v0i and Pi = P0i into the x and y components of 

equation (11c) gives 


Pk-t = i 


PKy = i 


(r^ + 4/3) kIV^ + (r^ + 1/3) n^KyVy + 

(r^ + 4/3) nlVy + (r^ + 1/3) n^KyV^ + nlVy 


\Ky^ 

AKj;, 


(B.4a) 

(B.4b) 


which can be combined using to obtain 


P = i (r^ + 4/3) [kxVx + nyVy 
= (r^ + 4/3) Cl/01 — Asin26' 
By substituting this into equation (11b) we obtain 

1 


2^-^A 




(B.5) 


Ci = 


+ 4/3 
P 

+ 4/3 


V ■ (iPAK0] 


(k, Ak) 01 


= — (k, Ak) Cl + 


A sin 26 


(k, Ak)0i, 


(B.6) 


+ 4/3 

where the inner product notation means (k, Ak) = This provides an expression 

* i 

for Cl/01 that, when substituted into (B.3), gives the growth rate 

(k. Ah) 


s = A(r(;'+ 4/3) (1 —0o)sin26* 


1 + (k. Ah) 


(B.7) 
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To understand the role of perturbation wavelength in controlling the growth rate we consider 
the anisotropy factor 

^ ^ {k,Ak) ^ [k\\ - (/C|| -/cx)sin^(^ + 0)] ^ 

1 + {k, Ak) 1 + [/c|I — (/c|| — /c_i_) sin^(6' + 0)] 

A rough measure of the effectiveness of anisotropy of permeability in controlling the angular 
dependence of the growth rate is the range of variation of F. The maximum and minimum 
values of F over all (0 + 0) are 


K^k\\ 

^ 1 +A’ 


_ 

““ 1 + K‘^ki_ ' 


We seek the stationary point of Fmax — TAin with respect to wavenumber k, 

^ d ^ ^ ^ ^ k\\ k± 

{l + k\\K^) {l + k^K^f 


= 


Taking a second derivative gives 

d^ 


{Fr, 


F ■ ) = 

2\2 max miny 


\/hkA 


2k‘^, 


2 /cm 


- . {l + k^K^ {l + k\\K^) 

which, evaluated at = 1/ ^/k\\k±, gives 


3 ’ 


^(k^)'3 (v+-'/h)<0. 

(yil+ v+) '' ’ 


(B.9) 


(B.IO) 

(B.ll) 


(B.12) 


(B.13) 


This result shows that this stationary point is a maximum in the total variation of F. Using the 
example of k\\ = 2, k± = 0.005 as in the numerical simulations, we hnd that the perturbation 
most affected by anisotropy has non-dimensional wavelength 2%/^ = 27r/v^T0 = 1.987... ~ 1, 
as expected. 
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